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ABSTRACT 

The problem of preconditioning the pseudospectral Chebyshev approximation 
of an elliptic operator Is considered* The numerical sensitiveness to 
variations of the coefficients of the operator are Investigated for two 
classes of preconditioning matrices: one arising from finite differences, the 
other from finite elements* The preconditioned system Is solved by a 
conjugate gradient type method, and by a DuFort-Frankel method with dynamical 
parameters* The methods are compared on some test problems with the 
Richardson method [12] and with the minimal residual Richardson method [17]* 
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Introduction 


Spectral approximations of elliptic boundary value problems lead to full 
and very ill-conditioned matrices* In the special case of constant 
coefficient operators, efficient direct methods have been proposed, [8], [91* 
For nonconstant coefficient problems, considerable attention has been devoted 
after Orszag's paper [12] to the simultaneous use of iterative methods and 
preconditioning techniques . 

In the present paper, we present and discuss the results of a number of 
numerical tests on the iterative solution of preconditioned systems arising 
from Chebyshev approximations. The first part is devoted to the analysis of 
the preconditioning of spectral matrices. The sensitiveness to variations of 
the coefficients, to leading and lower order terras is investigated. Besides 
the standard finite difference matrix proposed in [12] , we consider a finite 
element matrix, which essentially retains the same preconditioning properties, 
being moreover symmetric. In both cases, an incomplete factorization of 
Wong's type [16] is used* 

Two Iterative methods are considered next. A preconditioned conjugate 
gradient method (which has been recently used in fluid dynamics and transonic 
flow calculations via finite elements, see [5] and the references therein) 
resulted to be rather slow on the tested problems, although it may be very 
robust in more complicated situations. The DuFort-Frankel method (first 
applied by Gottlieb et. al. [6] , [7] to spectral calculations and here 

considered as a two-parameter preconditioned iterative method) yields good 
results when the optimal parameters are used. In order to overcome the 
difficulty of finding such parameters, we propose a modified version of the 
DuFort-Frankel method, devised according to a **mlnimal residual" strategy. 
The new method, compared with other iterative techniques in the literature, 
was the fastest in terms of speed of convergence. 
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No attempt Is made In this report to give theoretical justifications to 
the methods, nor to consider nonlinear problems* Both the aspects are, 
however, under investigation* 

Part of this work has been made while the authors were visiting ICASE* 
The numerical results reported here were obtained on the Honeywell 6040 at the 
the University of Pavla* Programs were written in double precision* The 
eigenvalues of Section 2 were obtained by EISPACK routines* 


2* The Preconditioning of ^ectral Matrices 

Let L be a smooth second-order elliptic partial differential operator 

1 2 2 
over the Interval ^ = (-1,1) or the square ^ =* (-1,1) * We consider 

homogeneous Dirichlet boundary conditions for L, i.e*, functions on which 

L acts will be assumed to vanish identically on the boundary* L will 

denote the Chebyshev pseudospectral approximation of L of order N* This 

means that the approximate solution is a polynomial of degree N and 

derivatives are computed after interpolating the function by a polynomial of 

degree N at the Chebyshev nodes ({x^ = cos > j**0,***,N if 

Q s q} • {x^,Xj}, 0 < l,j < N if Q ^ fi^)« We identify L^p with a matrix 

which maps the set of values of a polynomial u at the interior Chebyshev 

nodes into the set of values of the spectral approximation of Lu at the same 

nodes* 

It is known that L has a full structure. Moreover, its condition 

sp 

number is O(N^). These are considered negative aspects of spectral Chebyshev 
approximations versus finite difference and finite element methods* However, 
a tremendous improvement in the computational efficiency of spectral methods 
comes from the observation that L^p can be easily approximated by a sparse 
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matrix A, such that the condition number kTa ^ L 1 is close to 1 (see 

sp^ 

Orszag [12]), Throughout the paper, we refer to the ratio 

K = <(M) = |X |/|X ^ I as to the “condition number" of the matrix M, even 

max min 

when M is not symmetric* 

In the following, A will denote any matrix having these properties, and 
it will be called a preconditioning matrix. A is assumed to be related to 
some discretization of the operator L, usually by finite differences or 
finite elements • Sometimes we shall relate A to some other elliptic 
operator c^, with the same principal part as L. 

In one space dimension, the simplest way of building a preconditioning 
matrix is to use non-equally spaced finite differences at the Chebyshev nodes* 
The resulting matrix is tridiagonal, and it can be factorized in 0(N) 
operations* If Lu = "U^x* corresponding preconditioning matrix is given 

by A = where 

^ 2 ? 

"j - "j - "J+1 • 

In Table 2.1, the operator Lu = considered* The smallest and the 

largest eigenvalue X and X and the ratio k =» X /X ^ are 
° ° min max max min 

reported for both the matrices L^p and A"^ ^sp* 


( 2 . 1 ) 


jj hjhj.j 


-2 
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Table 2»1. Lu = 

Au = finite differences at Chebyshev points for Lu« 


N 

^sp 

A“*^ L 

sp 


X , 

X 

K 

A 

X 

K 


min 

max 


rain 

max 


4 

2.46 

.20 E2 

.80 El 

1. 

1.75 

1.75 

8 

2.47 

.21 E3 

.87 E2 

1. 

2.13 

2.13 

16 

2.47 

.32 E4 

.13 E4 

1. 

2.30 

2.30 

32 

2.47 

.50 E5 

.20 E5 

1. 

2.38 

2.38 

64 

2.47 

.80 E6 

.32 E6 

1. 

2.43 

2.43 

128 

2.47 

.13 E8 

.52 E7 

1. 

2.45 

2.45 


As expected, the largest eigenvalue of L grows like N^, while the 

sp 

eigenvalues of the preconditioned matrix A"^ L-_ lie in the Interval 

fap 

[!•, 2.5]. The spectrum of A"^ L exhibits a similar behavior even if the 

£>p 

elliptic operator L contains lower order terms (see [12]). 

The preconditioning properties of the matrix A seem to be rather 
insensitive to the lower order terms of L. The following table shows that 
the condition number k(A ^ ^sp^ kept small when A is just the finite 

difference approximation of the second-order terra of L, and the lower order 
terms are not prevailing* This implies that the preconditioning matrix can be 
kept fixed in solving nonlinear problems in which the lower order terms only 
change during the iterations* In all the cases considered below, the smallest 
eigenvalue is close to 1, and it converges to 1 from above as N 


increases 
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Table 2«2. Condition number k(A ^ L ) 

sp 

Lu = -u + 5u + Yu 

XX X 

Au = finite differences for = -u • 

XX 


N 

6=0. 
Y = 1. 

mi 

H 

■■ 

mm 

4 

1.38 

2.25 

— 

1.45 

2.29 

1.61 

8 

1.87 

3.45 

1.90 

2.01 

1.96 

16 

2.16 

4.30 

2.15 

2.28 

2.41 

32 

2.32 

4.73 

2.31 

2.47 

2.82 

64 

2.40 

4.92 

2.38 

2.85 

3.25 

128 

2.44 

4.99 

2.43 

3.08 

3.54 


When the magnitude of the lower order terms Is exceedingly large, the 
condition number of A"^ Lgp deteriorates. However the spectrum is still 
uniformly bounded in N. 


Table 2. 3. Condition number k(A ^ L 1 

sp" 

Lu = -u + 6u + Yu 

XX X 

Au = finite differences for = -u • 

XX 


N 

BHBI 

6 = 100. 
Y = 0. 


4 

3.95 

21.76 

217.43 

8 

17.70 

19.55 

195.39 

16 

28.90 

18.08 

181.05 

32 

35.89 

19.18 

177.83 

64 

39.17 

21.55 

176.24 

128 

40.57 

24.32 

204.05 
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A family of variable Coefficient operators Lu = with 

0 < < a(x) < , is considered in the next table. The eigenvalues of the 

matrix Lg^ are bounded independently of N, although the bound is larger 

than for the constant coefficient operator. The condition number < is close 
to the one in Table 2.1 when a moderate perturbation is applied, otherwise it 
grows slowly and linearly with the total variation of a. 


Table 2.4. Lu = - f(l + lo\^)u ) 

^ X^X 

Au = finite differences for Lu. 


N 

V = 0 

V = 1 

V = 2 


X K 

min 

X , < 

min 

X ^ K 

min 

4 

1.04 

2.49 

1.11 

5.03 

1.13 

7.09 

8 

1.01 

3.04 

1.10 

7.11 

1.01 

14.77 

16 

1.00 

3.27 

1.00 

7.70 

1.00 

21.65 

32 

1.00 

3.38 

1.00 

7.94 

1.00 

24.03 

64 

1.00 

3.43 

1.00 

8.06 

1.00 

24.44 

128 

1.00 

3.46 

1.00 

8.12 

1.00 

24.61 


VThen the coefficient a depends itself on the solution u (as in the 
full potential equation) one would not change the preconditioning matrix at 
each iteration, in order to save factorization time. This situation is 
simulated to a certain extent in the next table. The effects of 

preconditioning the spectral matrix of a variable coefficient operator by a 
constant coefficient operator matrix are reported. 















7 


Table 2.5. Lu = -((1 + 10^x^)u ) 

x^x 

Au = finite differences for S^\i » -u • 

XX 


N 

V = 0 

V =* 1 

^mln 

K 

^mln 

< 

4 

1.74 

2.23 

8.51 

3.70 

8 

1.48 

3.07 

4.70 

7.61 

16 

1.27 

3.79 

2.71 

12.79 

32 

1.16 

4.27 

1.83 

18.46 

64 

1.10 

4.60 

1.43 

23.29 


The spectrum of L__ Is bounded independently of N* < Is comparable 

with the one of Table 2.4 when the perturbation is moderate, but it becomes 
noticeably worse when the distance between the preconditioning and the 
spectral operators increases. In this case, if the factorization is carried 
out in 0(N) operations only (in one dimension or with an incomplete 
factorization) , the worsening of the condition number may not be balanced by 
the saving in factorization time (unless the computation of the entries of 
A is particularly expensive). 

The matrices A we considered so far arose from a finite difference 
approximation of the operator at the Chebyshev points. Even if ^ is 

formally self-adjoint, A is not symmetric, as well as Actually A 

splits up as A = D*A, with D diagonal and A symmetric. Some iterative 
techniques require the symmetry of the preconditioning matrix (see, e.g., the 
next section) . This can be accomplished by discretizing a suitable 

variational formulation of the elliptic operator via finite elements as 


follows 
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If ^\i = -(otu ) , the bilinear form associated to ^ is 


X X 


( 2 . 2 ) 


a(u,v) = / ot u^(vw)^ dx 


2 - I/9 

where w(x) = (1 - x ) The form a(*,*) is continuous and coercive on 

the weighted Sobolev space (cf* [2]), but it is not symmetric* 

However the “reduced” form 


(2.3) 


i(u,v) = Jt a u V w dx 


is still coercive and continuous on trivially symmetric- 
Assuming u and v continuous and piecewise linear between contiguous 
Chebyshev knots, we associate a matrix A = (2-3) by setting 


(2.4) 




where <|)^ Is continuous piecewise linear and ^k£* Instance, 

if ^u = -u , we have after dropping the common factor ^ 


N • 


(2.5) 


jj 


2 2 


j.j-1 


-1 


j-1 


j,j+l 




(compare with (2.1)). The spectrum of A behaves like the spectrum of the 
corresponding finite difference matrix, and the preconditioning properties are 
only slightly worse, as shown in the next table. 
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Table 2% 6. Lu = "^xx 

Au =» finite elements at the Chebyshev 
points for Lu (see(2.5))* 


N 

X (a"^ l ) 

mln^ sp"^ 

X (a“^ l ) 

max ^ sp 

k(a“^ L ) 
^ sp^ 

4 

1.25 

3.29 

2.63 

8 

1.16 

4.10 

3.55 

16 

1.13 

4.53 

3.99 

32 

1.13 

4.74 

4.19 

64 

1.13 

4.84 

4.29 

128 

1.13 

4.89 

4.33 


Up to now we considered one-dimensional problems. In 2D one can still use a 
finite difference or finite element matrix, say B, in the preconditioning. 
The corresponding results are similar to those in ID. However, the exact 
'’inversion” of such a matrix is more expensive, since the factors in its LU 
decomposition have a bandwidth of order 0(N) instead of 0(1). In order to 
overcome this drawback, different techniques of incomplete factorization have 
been successfully proposed (cf., e.g., [10], [11], [16]). The idea is to 
replace the exact factors L and U by some approximations L and U of 
them, which retain a very sparse structure. L and U are computed by 

incomplete steps of Gaussian elimination, under the condition that certain 
quantities depending on the product L U agree with the corresponding 
quantities for B. The matrix A « L U is then used in the preconditioning. 

In our computations, the inqoraplete factorization was done according to 
Wong^s row-sums agreement condition ([16]). Namely, let b^®^ and b^^^ 
denote the diagonal and the off-diagonals of a m-th order matrix B, l.e.. 
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= {b^ 1+k^^ ^ i,t+k < m} . If B is the finite difference matrix for a 
second-order operator at the Chebyshev points in the square, then only b^^\ 
b^^^ ^ , and are not Ident ically zero. The Incomplete factors L 

and U of B have and u^^\ respectively 

as nonzero (of f )-diagonals* is chosen to be ^ while the off- 

diagonal elements are easily determined by the condition that a^^^ ^ = b^^^^ 
and a^^^ ^ b^^\ where A “ L IJ. Finally, is such that the sum of 

each row in A equals the corresponding sum in B, 

Henceforth, we list some results about the preconditioning by an 
Incomplete factorized finite difference matrix (for other results see [17]). 
Table 2.7 refers to constant coefficients operators. The different ratios 
between the coefficients of u„^ and u,_„ are supposed to mimic the effect 

AX yy 

of the Stretching of coordinates in a mapping process. The spectral matrix of 
a variable coefficients operator was preconditioned by the finite differences 
representation of the same operator (Table 2.8), or by that of a constant 
coefficient operator (Table 2.9). 


Table 2. 7. Lu = -uu - u 

XX yy 

Au - incomplete factorized finite 
difference matrix for LU. 


N 

a = 1. 

a = 10. 

a = 100. 


A ^ K 

min 

X K 

min 

X K 

rain 

4 

1.08 

1.72 

1.01 

1.75 

1.00 

1.76 

8 

1.06 

2.72 

1.03 

2.43 

1.01 

2.13 

16 

1.04 

4.06 

1.03 

5.34 

1.01 

3.27 
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Table 2. 8. 

Lu=(otu) - (Pu) 

^ x^x y y 

a = 3 = 1 + lOx^y^ 

Au ** incomplete factorized finite 
difference matrix for Lu. 


Table 2.9. 

Lu as in Table 2.8 
Au = incomplete factorized 
finite differences for 
^\x = - u. 


XX 


yy 


N 


K 

min 


4 

1.48 

6.89 

8 

1.44 

10.94 

16 

1.23 

19.90 


N 

X 

K 

min 


4 

1.09 

3.29 

8 

1.08 

4.92 

16 

1.04 

9.33 


Unlike the case of complete factorization, the condition number grows linearly 
with the number of unknowns. However, it ranges within moderate bounds 
(except when a different operator is used in the preconditioning). This gives 
evidence to the convenience of using incompletely factorized preconditioning 
matrices in spectral calculations. Better results can be achieved, with 
slightly more computational effort, by a higher order incomplete factorization 
in which L and IJ have one more nonzero off-diagonal (see [16] for more 
details) • 


3. A Conjugate Gradient Method 

Even if the differential operator L is self-adjoint, the matrix arising 
from a Chebyshev spectral approximation is not symmetric. Thus, one can apply 
the standard conjugate gradient method (CG) to the normal equations of the 
preconditioned system; or one can use CG-type methods for nonsymmetrlc 
systems, like those proposed by Vinsome [13], Young and Jea [14], Axelsson 
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[1], or those by Concus and Golub [3], Wldlund [15]. The methods of the first 
class may require the storage of back steps of the solution (see however, Wong 
[17] for an application of a truncated version of [1 ] to spectral 
calculations), while the methods of the second class need that the symmetric 
part of the system be easily solvable. 

In the previous section it has been pointed out that the spectral matrix 
can be preconditioned using a symmetric positive definite matrix, connected 
with some finite element approximation of the elliptic operator. This 
suggests the use of the following preconditioned version of the CG method 
(see, e.g., [5]): Minimize 


(3.1) 


J(u) = r^ r 


r = L u - f 
sp 


by CG iterations in I?' equipped with the inner product 
(3.2) ((u,v)) = u*^ Av. 

The corresponding algorithm is as follows. 


Given vP c compute z^ = A ^(f - u^) 

0 -1 ^T 0 0 0 

g = A z , w =* g - 


Then set for k = 0, !,•••: 


(3.3) / 


k+1 

k , k 

k 

u 

= u + ot 

w , 

k+1 

k k 

*-l 

z 

= z - ot 

A 

k+1 

tT 

k+1 

g 

= A L 

sp 

z 

k+1 

k+1 , . 

^,k+l 

w 

=* g + 

Y 


where ot = 




[l w^,A ^ L w^) 
^ sp sp 


sp 


k+1 k k+1 . 
^ - R »R 

((g’",g^)) 
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T ^ T k+1 

Here (u,v) = u v denotes the Euclidean product in Ik** The product L z 

sp 

can be executed through Fast Fourier Transforms, and the entries of the matrix 
Lgp need not to be computed. Actually, assume that 
the Chebyshev pseudospectral approximation of Lu » 

the N-th degree polynomial interpolating w at the nodes x^ , j«0, • • • ,N. We 

identify a N-th degree polynomial vanishing at x = ±1 with the vector of its 
values at Xj , j=l, • • • ,N-1. Recall that 

1 N-1 

(3.4) / u(x) v(x) W(x) dx = TT I u(x ) v(x ) + {u(-D v(-l) + u(l) v(l)} 

-1 W J J N 


for any u,v such that uv e 1^2N-1* Then 

N ^ 

(u,L^ v) = (l u,v) ~ / [P )] V w dx. 

sp ^ ^ sp ’ 1 m_i ^ c X ‘‘x 

^ W ^ -1 


Integration by parts and several applications of (3.4) yield 


(3.5) 


V = -fp 2 ) + 32 

sp ^ C X 

where z = ot(v + Bv) , 

X 


(0 

6(x) “ •;;p (x) = — - — 2 
1 - X 


Similar expressions hold in two dimensions « 


Algorithm (3.3) was used to compute the spectral solution for the test 
problems : 

Lu = - (otu 1 =* f, -1 < X < 1 

^ x^x / 

a = 1. or ot(x) = 1 + lOx^ 
u(x) = sin TTx, 
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and 


(3.7) 


/ Lu = - [otu 1 - (au ) = f -1 < x,y < 1 

I ^ x^x 

I 2 2 

/ a = 1, or ot(x,y) = 1 + lOx y 

\u(x,y) = sin TTx sin iry. 


In the following tables, we report the minimum number NIT of Iterations 
required to get RES < 10 , where the relative residual Is defined by 


(3.8) 


RES 


2 - (r.r) 

(f,f) ’ 


X = 



u. 


The Initial guess was u^ = 0. ERR Is the corresponding relative error on 
the solution 


(3.9) 


ERR = 


II u - u 

sp exact 

ilu ^11 

exact 


\u 2 

where Wuil = (u,u) ^ Is the discrete ^ -norm on the grid. 


Table 3.1. CG Method for Problem (3.6) 

Au = finite element matrix for Lu. 


N 

a = 1. 

0=1+ lOx^ 


NIT 

RES 

ERR 

NIT 

RES 

ERR 

4 

1 

.31 E-17 

.18 EO 

1 

.90 E-17 

.11 E-l 

8 

3 

.12 E-16 

.31 E-3 

3 

.19 E-13 

.40 E-3 

16 


.53 E-14 

.27 E-11 

8 

.13 E-15 

.61 E-11 

32 


.58 E-8 

.12 E-9 

16 

.60 E-12 

.88 E-14 

64 

20 

.52 E-8 

.84 E-10 

24 

.24 E-8 

.18 E-10 

128 

26 

.52 E-8 

.47 E-11 

29 

.35 E-8 

.13 E-10 
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Table 3. 2* CG Method for Problem (3.7) 

Au = incomplete factorized finite element matrix for Lu. 


N 

P 

ill 

. 

a 

« 1 + lOx^y^ 

1 


NIT 

RES 

ERR 

NIT 

RES 

ERR 

4 

21 

.31 E-8 

.18 EO 

32 

.88 E-8 

.10 E-1 

8 

44 

.75 E-8 

.31 E-3 

72 

.99 E-8 

.15 E-3 

16 

80 

.99 E-8 

.89 E-9 

70 

.13 E-2 

.56 E-4 


It is seen that the number of iterations NIT to match the stopping criterion 
RES < 10 increases sub linearly in ID and linearly in 2D with the degree N 
of polynomials. This seems qualitatively in accordance with the behavior of 
the condition number of the matrix ^sp* however, we are unable to find a 
satisfactory explanation to the slow convergence properties of the method in 
2D. 


4. The DuFort-Frankel (DF) Method 

The DuFort-Frankel method can be applied to the numerical solution of 
steady-state equations 

(4.1) Bu = g, 

(the eigenvalues of B having positive real parts) as an iterative procedure 
depending on two parameters 6 and y: 
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(4.2) 


k+1 k-1 

u - u 

26 


„ K k+1 ^ K , K-1^ 

Bu - y[u - 2u + u J. 


k-1 


This can be written as a one-step method in the form 


k+11 


r k 1 


- 

u 

k 

= G(B;6,y) 

u 

k-1 

1 + 26y 
26 

g 

_u 


_ u 


0 


with proper definition of the matrix G. 

The DF scheme was studied in connection with spectral methods by Gottlieb 
and Gustafsson [6] and by Funaro [4]. If B has real strictly positive 
eigenvalues (the largest and the smallest eigenvalues being denoted by 
and respectively) , then it is seen that the method is convergent if 


(4.4) 


Y > Y. 


LIM 


X 

max 

4 • 


Moreover, Funaro [4] proves that the spectral radius P(G) as a function of 
6 and Y has a curve of local minima (with respect to increments in the 6 
or in the Y direction separately) given by the branches of hyperbola 


(4.5) 


1 + 


Y = 


max 




max 


If 


X ^ + X 

min max 

4 


(4.6) 


Y = 


1 + 6^X^ 

min 


46^X 


min 


If 


Y > 


X , + X 

min max 

4 


p(G) attains Its absolute minimum at the intersection of the two branches, 
l.e., at the "optimal parameters" 


(4.7) 


6* = 


* _ min max 
Y * A 


VX~X 
min max 
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where 


(4.8) 


P(G) 


opt 


* 

P 


8 


r 


max 


r 


max 


/A 

A 


min - 1 


min + 1 


The DF method with the optimal parameters (4.7) was applied to the 
solution of the test problems (3.6) - (3.7) by a preconditioned spectral 

method. Hence, we set in (4.2) Bu = ^sp^ 8 “ where A is 

the finite difference matrix associated to L, incompletely factorized in 

2D. One DF iteration requires one multiplication z =* Lgp u^ and one 
forward-backward substitution Aw =* z. The optimal parameters were computed 
using the exact values of A . and A obtained in the previous 

section. The initial guess was u^ = 0, while u^ was computed by a step of 
the Modified Euler method for the preconditioned system. NIT, RES, ERR are 
defined as in Section 3, formulae (3.8) - (3.9). 


Table 4.I. DF Method with optimal parameters for Problem (3.6) 
Au « finite differences for Lu. 


N 

a = 1. 

a = 1 + lOx^ 


NIT 

RES 

ERR 

NIT 

RES 

ERR 

4 

9 

.33 E-8 

.18 EO 

19 

.41 E-8 

.11 E-l 

8 

12 

.19 E-8 

.13 E-3 

24 

.97 E-8 

.40 E-3 

16 

12 

.98 E-8 

.39 E-8 

26 

.52 E-8 

.56 E-8 

32 

14 

.25 E-8 

.14 E-8 

29 

.37 E-8 

.79 E-8 

64 

14 

.29 E-8 

.29 E-8 

28 

.83 E-8 

.75 E-8 

128 

14 

.57 E-8 

.56 E-8 

28 

1 

.88 E-8 

.86 E-8 
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Table 4»2. DF Method with optimal parameters for Problem (3*7) 

Au « incomplete factorized finite differences for Lu* 


N 

a = 1. 

a 

= 1 + lOx^y 

2 


NIT 

RES 

ERR 

NIT 

RES 

ERR 

4 

9 

.70 E-8 

.18 EO 

15 

.20 E-8 

.10 E-l 

8 

14 

.97 E-8 

.13 E-3 

20 

.84 E-8 

.15 E-3 

16 

20 

.56 E-8 

.11 E-8 

30 

.98 E-8 

.64 E-8 

32 

59 

.86 E-8 

.48 E-8 

49 

.82 E-8 

.92 E-9 


It is seen that the number of iterations needed to satisfy the stopping test 
is bounded as a function of N in the ID tests, while it is linearly growing 
in 2D, This corresponds to the behavior of the condition number of the matrix 
reported in Tables 2.4 and 2.8. 

sp 

Moreover, NIT is comparable with the one relative to the CG method in ID, 
and definitely smaller in 2D. Since one DF iteration is faster than one CG 
iteration (by a factor of 1.7 both in ID and in 2D) we conclude that the DF 
method with optimal parameters exhibits a globally better behavior than the CG 
method on the tested problems. 

The speedup in the convergence due to the use of a preconditioning 

technique is particularly impressive for the DF method. This is suggested by 

formula (4*8), which shows the dependence of the optimal spectral radius on 

the condition number of B. Table 4.3 reports the performance of the DF 

method with optimal parameters without preconditioning (i.e., Bu =» u) for 

sp 

problem (3.6) with a = 1. (compare with Table 4.1) 
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Table 4.3. DuFort-Frankel method without preconditioning 


N 

4 

8 

16 

32 

64 

NIT 

23 

84 

327 

>400 

»400 


The practical Interest of formulae (4.7) relies on the explicit knowledge 

of X and X , which is rarely the case* Approximate values of 5 
min max' ^ 

* 

and Y may be obtained in different ways, for instance by estimates on the 

eigenvalues of B or by extrapolation of correct values of 6* and Y* 

computed on coarser grids. It was found that linear extrapolation on the 

* 

parameters as functions of N may lead to negative values of Y • Instead, 
linear extrapolation on the ratios of contiguous values of the parameters 
gives accurate answers. Actually, the case N = 32 in Table 4.2 was run with 
extrapolated "optimal** parameters. 

Unfortunately, the method appears to be rather sensitive to the choice of 

parameters, especially around the curve of optimality. The qualitative 

behavior of P(G) as a function of Y for fixed 6 (or conversely) is 

similar to the one encountered in a SOR method. Table 4.4 shows the values of 

2 2 

NIT for problem (3.6) with a 1 + lOx y , N = 32 and finite differences 

preconditioning, as a function of Y and 5. 
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Table 4.4. NIT as function of Y and 6* 


‘ x6 

4 

»400 

193 

438 

»400 

»400 

2 

>400 

98 

162 

416 

»400 

1 

>400 

49 

97 

183 

392 

V2 

>400 

87 

61 

134 

196 

V4 

»400 

147 

145 

86 

151 


Y /y* 

lim' ' 

1 

2 

4 

8 

► 


The previous considerations suggest that the DF method, although in 
principle very powerful, may be poorly efficient in applications if the user 
attempts to fix the constants 5 and Y once for all. 

However, it is possible to transform the DF method into a completely 
parameter-free iterative scheme following a "minimal residual" strategy which 
has been proven successful in connection with other Iterative schemes. The 
parameters Y and 6 are computed at each iteration in order to minimize the 
£^-norra of the residual r = f - L_„ u. Given u^,r^ and u^“^, then 

sp 

and r^"^^ are defined according to (4.2) 

(4.9) 
where 


k+1 .-1 k , k , k-1 

u = Cj^ A r + C2 u + Cg u 


k+1 , k . k , k-1 

r “ -c, L A r +c^r+c^r 
1 sp 2 3 


25 ^ 46y 

" 1 + 26y ’ “ 1 + 26y 


(4.10) 


» 


1 - 25y 
“ 1 + 26y * 
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is minimized if one sets in (4.10) 


(4.11) 


where p =* q =« Lgp r^, s = + p, a » (p,r^“^) /(p,s) and 

P * (p,q)/(p,s) . 

This algorithm can be called ’’Minimal Residual DuFort-Frankel” (MRDF) 
Method. One MRDF iteration requires one forwarcf^backward substitution 
Az = r^ and one multiplication w = z; moreover, needs to be stored 

with u^"“^. Note that if u^ = u^ but r^ ^ 0 the algorithm cannot 

converge. Hence u^ should be chosen in such a way that u^ - u® and r® 
be roughly comparable. For Instance u^ can be computed from u^ by one 
step of the Minimal Residual Richardson method (see Section 5 (b)). 

Tables 4.5 and 4.6 are analogous of Tables 4.1 and 4.2, except that the 
MRDF method was used instead of the DF method with optimal parameters. 

In ID, the gain in the speed of convergence over the DF method with optimal 
parameters is spectacular, although this may depend on particular 

circumstances. In 2D, the improvement of the performances is less impressive. 

However, one must not forget that the main Improvement of the MRDF over 
the DF method relies on the complete automatization in the choice of 
parameters. 


5*^ = i. - ”s) 

2 (q.q - 3s) 


1 (p»26^q - r^~^) 


26 ' 


(p.s) 
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Table 4«5. MRDF Method for Problem (3.6) 

Au = finite differences for Lu. 


N 

a = 1. 


a = 1 + lOx 

2 


NIT 

RES 

ERR 

NIT 

RES 

ERR 

4 

1 

.11 E-17 

.18 EO 

1 

.70 E-17 

.11 E-l 

8 

5 

.48 E-9 

.13 E-3 

8 

.55 E-8 

.40 E-3 

16 

7 

.71 E-8 

.38 E-9 

11 

.46 E-8 

.35 E-8 

32 

4 

.56 E-9 

.13 E-9 

9 

.48 E-8 

.22 E-8 

64 

3 

.10 E-9 

.20 E-10 

3 

.45 E-8 

.18 E-8 

128 

2 

.33 E-9 

.68 E-10 

2 

.10 E-8 

.56 E-9 


Table 4«6. MRDF Method for Problem (3.7) 

Au = incomplete factorized finite difference matrix for Lu. 


N 

a = 1. 

a = 1 + lOx^y 

2 


NIT 

RES 

ERR 

NIT 

RES 

ERR 

4 

7 

.23 E-8 

.18 EO 

10 

.84 E-8 

.10 E-l 

8 

13 

.47 E-8 

.13 E-3 

20 

.77 E-8 

.15 E-3 

16 

19 

.63 E-8 

.80 E-9 

29 

.90 E-8 

.79 E-8 

32 

36 

.50 E-8 

.96 E-9 

46 

.78 E-8 

.86 E-9 


5* Comparisons with Other Methods 

The preconditioned CG and DF methods were compared with two other 
iterative techniques recently suggested for spectral calculations: the 

Richardson Iteration proposed by Orszag [12], and the Minimal Residual 






























Richardson method proposed by Wong [17]. We briefly review these techniques 
and we report for the sake of completeness their behavior on the test problems 
used throughout this paper. 
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(a) Richardson Method ([12]) : 

Given u®, compute u^"*"^ from u^ by solving 

(5.1) = Au^ - oi(l u^ - f), 

sp 

where 0 < a < 2/X , X and X being the smallest and the largest 

max min max ” ° 

eigenvalue of A"^ The optimal value of ot, 

sp 


(5.2) 


a 

opt 


2 


X ^ + X 

min max 


was computed exactly and used in the following tests. One iteration requires 
one multiplication z = w and one forward-backward substitution Ax = b. 


Table 5.I. Richardson Method for Problem (3.6) 
Au = finite differences for Lu. 


N 

a = 1. 


0 = 1 + 10x‘ 



NIT 

RES 

ERR 

NIT 

RES 

ERR 

4 

8 

.90 E-8 

.18 EO 

33 

.63 E-8 

.11 E-l 

8 

17 

.81 E-8 

.13 E-3 

62 

.95 E-8 

.40 E-3 

16 

20 

.77 E-8 

.54 E-8 

71 

.82 E-8 

.68 E-8 

32 

21 

.70 E-8 

.64 E-8 

73 

.96 E-8 

.88 E-8 

64 

22 

.70 E-8 

.41 E-8 

74 

.97 E-8 

.94 E-8 

128 

22 

.42 E-8 

.51 E-8 

75 

.87 E-8 

.86 E-8 
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Table 5.2. Richardson Method for Problem (3# 7) 

Au = incomplete factorized finite difference matrix for Lu. 


N 

a = 1. 

a = 1 + lOx^ 

2 

y 


NIT 

RES 

ERR 

NIT 

RES 

ERR 

4 

12 

.87 E-8 

.18 EO 

23 

.84 E-8 

.10 E-1 

8 

24 

.71 E-8 

.13 E-3 

45 

.69 E-8 

.15 E-3 

16 

39 

.92 E-8 

.92 E-8 

90 

.99 E-8 

.21 E-8 


(b) Minimal Residual Richardson (MRR) Method ([17]); 

In the previous scheme, compute a = a at each iteration in order to 
minimize the residual Hence one gets: 


Given u , compute r = f - u , z = A r , 


then set 

k+1 k , k k , k 

u = u + cx z where a = 


(r ,L 




(5-3) 


(l z ,L z ) 
^ sp sp ^ 




k+1 k k - k 
r = r - ot L z 
sp 


k+1 ,-l k+1 

z = A r 


The computational effort per iteration is comparable with that of the 
Richardson method. Note that this method is obtained from the previous one by 
the same strategy used in deriving the MRDF from the pure DF method. 

















Table 5.3. MRR Method for Problem (3.6) 

Au = finite differences for Lu 


.12 E-17 
.32 E-8 
.78 E-8 
.56 E-9 
.14 E-9 
.34 E-9 


.18 EO 
.13 E-3 
.29 E-9 
.19 E-11 
.12 E-10 
.48 E-10 


a = 1 + lOx 


RES 

ERR 

.15 E-17 

.11 E-l 

.22 E-10 

.40 E-3 

.76 e-8 

.28 E-8 1 

.62 E-8 

.37 E-8 

.58 E-8 

.15 E-8 

.16 E-8 

.12 E-8 


Table 5.4. MRR Method for Problem (3.7) 

Au = incomplete factorized finite difference matrix for Lu 


1 + lOx y 


NIT 

RES 

ERR 

NIT 

RES 

ERR 

9 

.98 E-8 

.18 EO 

18 

.96 E-8 

.10 E-l 

18 

.11 E-8 

.13 E-3 

22 

.29 E-8 

.15 E-3 

23 

.90 E-8 

.53 E-8 

32 

.14 E-8 

.60 E-9 

58 

.88 E-8 

.19 E-8 

58 

.89 E-8 

.43 E-9 


(c) Comparisons 

The speed of convergence of the methods previously discussed was compared 
on the basis of the number of iterations and the CPU time. Two significant 
cases were considered. 
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CASE 1; Problem (3.6) with a = 1 + lOx^ 

N = 128, i.e*, 127 grid points in the interval (-1,1)* 

CASE 2: Problem (3.7) with a = 1 + lOx^y^ 

N = 32, l.e., 31 ^ 31 grid points In the square (-1,1)^. 

Define for the sake of simplicity the following labels: 

A Richardson method (5.1) with optimal parameter (5.2) 

B Minimal residual method (5.3) 

C Conjugate gradient method (3.3) 

D DuFort-Frankel method (4.2) with optimal parameters (4.7) 

E Minimal residual DuFort-Frankel method (4.9) 

We used the standard finite difference (finite element for method C) 

preconditioning matrix on the spectral grid, incompletely factorized in Case 2 

according to Wong's method described in Section 2. The optimal parameters 

were computed with the exact values of A and A , = 0 was the 

^ min max 

initial guess. 

The results in Figure 5.1 and in Figure 5.2 are in a sense machine- and 
programmer- independent. The relative performances of the methods can be 
analyzed according to the global CPU-tlrae consumption, using the following 
table. 


Table 5% 5. CPU-time per iteration in hours 


METHOD 

A 

B 

c 

D 

E 

CASE 1 

.272 E-3 

.280 E-3 

.467 E-3 

.273 E-3 

.285 E-3 

CASE 2 

.128 E-Z 

.128 E-2 

.217 E-2 

.127 E-2 

.130 E-2 














Hence, Figure 5.1 and Figure 5.2 also summarize the relative performances of 
the methods in terms of cost, except for method C which is roughly 1. 7 times 
slower than the others. 
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Coizmients 

Globally, the results confirm the utility of preconditioning techniques 
in spectral calculations : few iterations are needed to reach the spectral 

accuracy, which corresponds in the test problems to a relative residual of 

10 -> 8 . 

Methods A and D behave exactly like predicted by the theory; the error 

K—1 

is reduced at each Iteration by a factor — rr for method A, and ~ for 

/ST 

method B (k Is the condition number of the preconditioned matrix). 

The conjugate gradient method gives contradictory answers in terms of 
speed of convergence: in ID the factor of reduction of the error is smaller 

than that for method D, while in 2D it is comparable with that of method A. 
In both cases, the method turns out to be not competitive in terms of computer 
time. 

The "minimal residual" strategy is always winning over the "optimal 
parameter" strategy, also where the exact optimal parameters can be used. In 
particular, the MRR method is superior even to the Richardson method with 
Chebyshev acceleration, proposed in [12] (according to [12] , p. 86, the 
Chebyshev acceleration increases the speed of Richardson method by a factor of 
2, although it requires the extra-storage of the vector u^“^). 

The MRDF method requires the storage of u^*"^ and r^"^, being a two- 
step method. However, the extra memory required results in a better accuracy, 
and the MRDF method appears in all cases the fastest method among those tested 
in this report. 
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